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We consider parameter estimation for the spread of the Neolithic 
incipient farming across Europe using radiocarbon dates. We model 
the arrival time of farming at radiocarbon-dated, early Neolithic sites 
by a numerical solution to an advancing wavefront. We allow for 
(technical) uncertainty in the radiocarbon data, lack-of-fit of the de- 
terministic model and use a Gaussian process to smooth spatial devi- 
ations from the model. Inference for the parameters in the wavefront 
model is complicated by the computational cost required to produce 
a single numerical solution. We therefore employ Gaussian process 
emulators for the arrival time of the advancing wavefront at each 
radiocarbon-dated site. We validate our model using predictive sim- 
ulations. 

1. Introduction. The Neolithic (10,000-4500 BC), the last period of the 
Stone Age, was the dawn of a new era in the development of the human race. 
Its defining innovation, the transition from foraging to food production based 
on domesticated cereals and animals, was one of the most important steps 
made by humanity in moving toward modern societies. Among dramatic 
changes resulting from the advent of farming were sedentary living, rapid 
population growth, gradual emergence of urban conglomerates, division of 
labor and the development of complex social structures. Another important 
innovation of the period was pottery making. The mechanism of the spread 
of the Neolithic remains an important and fascinating question. The most 
striking large-scale feature of this process is its regular character, whereby 
farming technologies spread across Europe at a well-defined average speed of 
about U = 1 km/year. Impressive and convincing evidence for this emerged 
as soon as radiocarbon ( 14 C) dates for early Neolithic sites became available 
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Fig. 1. Left-hand panel: the location of 14 C-dated early Neolithic sites used in our anal- 
ysis and their age (calibrated BC, color coded). Right-hand panel: the dimensionless dif- 
fusivity in the wavefront model (5), with a mild linear gradient toward south, local 
topographic variations, and the gradual decline to zero away from the coastlines. 



[Ammerman and Cavalli-Sforza (1971)], and later studies have confirmed 
the early results [Gkiasta et al. (2003)]. Figure 1 displays the 1 C dates used 
in our analysis, and they clearly show the gradual, regular spread of the 
Neolithic over a time span of over 4000 years from its origin in the Near 
East around 9000 years ago, to the north-west of Europe. 

It is worth emphasising that "the Neolithic" is not a single phenomenon 
and, although it is often characterized by archaeologists as a "package" of 
related traits, the individual elements of the package need not have been 
transmitted simultaneously. This might suggest the need to model each dis- 
tinct element separately. However, several elements do appear to have trav- 
eled together during the spread of farming into most of Europe [Burger and 
Thomas (2011)] and the local transition to the full package may have been 
more rapid than has often been assumed [Rowley-Conwy (2011)]. Therefore, 
we choose to model the Neolithic spread as a single phenomenon. 

Despite the overall regular character of the expansion, there are notable 
regional variations in the spread of the Neolithic [Gkiasta et al. (2003), 
Bocquet-Appel et al. (2009)]. Careful analysis of the radiometric and archae- 
ological evidence has identified various major perturbations: First, slowing 
down at topographic obstacles such as mountain ranges. Second, latitudi- 
nal retardation above 54°N latitude [Ammerman and Cavalli-Sforza (1971)], 
presumably due to competition with the pre-existing Mesolithic population, 
whose density was especially high in the north [Isern and Fort (2010)]; the 
negative influence of the soil type and harsher climate in the north also can- 
not be excluded. Third, significant acceleration along the Danube-Rhine cor- 
ridor and the west Mediterranean [Ammerman and Cavalli-Sforza (1971)]; 
see also Davison et al. (2006). These variations in the rate of spread com- 
plicate quantitative interpretations of the 14 C data. 
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Moreover, the data themselves are far from being complete and accurate. 
There are inherent errors from the radiocarbon measurement in the labora- 
tory [Scott, Cook and Naysmith (2007)] and uncertainties in the conversion 
of the 14 C into the calendar age, known as the calibration of the 14 C dates 
[Buck et al. (1991), Blackwell and Buck (2008)]. Further errors result from 
archaeological factors, such as uncertain attribution of the dated object to 
the early farming activity, disturbed stratigraphy of archaeological sites, etc. 
[Aitken (1990)]. Perhaps even more importantly, objects related to the first 
appearance of farming at a given location are not necessarily related to the 
first arrival of farming to its larger region. The occupation of many Neolithic 
sites could be a secondary process behind the wave of advance. The only ob- 
vious method to identify such sites is by continuity with its neighbors. 

Given the obviously random nature of the demographic processes un- 
derlying the Neolithic expansion, and the abundance of both random and 
systematic errors in the quantitative (mainly radiometric) evidence avail- 
able, it is clear that serious statistical analysis is required to quantify the 
spread of the Neolithic. All previous work in this direction has determined 
the parameter selection for various models and assessed goodness of fit by 
using parameter scans. To the best of our knowledge, the present paper is 
the first to use formal statistical inference methods to determine appropriate 
ranges for model parameters and assess the fit of the model. 

We use a carefully selected set of C dates for the earliest Neolithic sites 
across Europe to fit a model of the Neolithic expansion. The model relies 
on the fundamental fact, established empirically as briefly explained above, 
that the incipient farming spread at a nearly constant average speed. We 
also include regional variations in the speed associated with the latitudinal 
gradient, topography and major waterways. Our statistical model adopts an 
error structure which accounts for both measurement error in the radiocar- 
bon dating process and a spatial error process which smooths departures 
from the wavefront model. We describe how to obtain the Bayesian poste- 
rior distribution for the unknown model parameters using simulations from 
the deterministic model. Simulations of the expanding Neolithic front are 
computationally expensive, even at a prescribed deterministic velocity, pre- 
cluding their use for Bayesian inference. Therefore, we construct a Gaussian 
process emulator [Santner, Williams and Notz (2003)] for the model, that is, 
a stochastic approximation to the arrival times obtained from this model. 
Such methods are widely used in the computer models literature; see, for 
example, Kennedy and O'Hagan (2001) and references therein. Rather than 
build a complex space-time emulator to describe the wavefront model, we 
pragmatically build individual emulators for each site and rely on the spatial 
Gaussian error process to smooth across geographical space. Finally, we fit 
this emulator model to our data and describe the inferences we draw. 
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The remainder of this paper is organized as follows. We describe the C 
data used in our analysis in Section 2. The implementation of the wavefront 
model and the statistical model which links the wavefront model to the 
data are presented in Section 3. The Gaussian process emulator model and 
its use in constructing an inference scheme are the subject of Sections 4 
and 5. Finally, the results of fitting our model to the radiocarbon data are 
reported in Section 6, and we conclude with a discussion in Section 7. 

2. Selection and analysis of radiocarbon dates. It is clear from the dis- 
cussion in Section 1 that 14 C data need to be carefully selected for model 
parameter inference. Databases of radiocarbon dates contain several thou- 
sand entries for Neolithic sites across Europe [e.g., Gkiasta et al. (2003)], 
but a relative minority of them are suitable for the studies of the initial 
spread of the Neolithic. In this paper, we use carefully selected 14 C dates 
from 302 earliest Neolithic sites in Southern and Western Europe from the 
compilation of Davison et al. (2009), where a detailed discussion of the se- 
lection criteria can be found. The objects dated include grains and seeds, 
pottery, bone, shells and animal horns, wood, charcoal and peat, and soil. 
We note that, in general, there may be issues in using such a heterogeneous 
data set, such as those relating to "chronometric hygiene" and the problems 
that may arise with certain materials or treatments; see, for example, Pettitt 
et al. (2003) and the discussion in Steele (2010). However, we are not aware 
of any issues with the provenance of the samples we use in our analysis; see 
Davison et al. (2009) for further details. 

The putative origin of the expansion is in an extended region in the Near 
East, and, following the results of Ammerman and Cavalli-Sforza (1971), 
later confirmed by many authors, it is reasonable to place it near Jericho. 
Similarly to Davison et al. (2009), the starting date and position are adopted 
as 6572 BC and (</>, A) = (41.PN, 37.1°E), the Tell Kashkashok site near 
Jericho. 

Many sites in the selection have several dates (often 3-10, and occasionally 
30-50) which can be treated as the arrival time contaminated by errors. Each 
date is published together with the accuracy of radiocarbon measurement in 
the laboratory, but we recall that other sources of random and systematic 
errors are present as well. Suppose that, at site i, we have rxii calibrated 
14 C dates (j = 1, . . . , rrn, i = 1, . . . , n) and their standard deviations aij. 
These data are calculated from the central 99.7% interval of the calibrated 
probability distributions for each object, obtained using calibration curve 
IntCal04: Reimer et al. (2004). Specifically, each tij is taken to be the mid- 
point of its calibrated interval and aij to be one-sixth of the length of the 
interval. This data set can be found in the supplementary material [Bagga- 
ley et al. (2013)] for this paper. We can calculate summaries for the data at 
each site, namely, an arrival time ij and its accuracy <7j, as defined in this 
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context by Scott, Cook and Naysmith (2007), by appropriately weighting 
the data, as 

mi mi mi 

(1) *i = I>i^ 2 /X>y 2 and o? = l/X>y 3 . 

j=i j=i j=i 

Clearly, there is some loss of information in not using the full calibrated 
probability distribution for each object. However, to include this within our 
analysis by using, for example, a normal mixture model representation for 
each distribution, would make our analysis considerably more complicated. 
In any case, the model we do use in our analysis is one for the weighted means 
ti, whose calibrated distributions will be roughly normally distributed due 
to the central limit theorem. The left-hand panel of Figure 1 shows the sites 
on the map, color-coded according to the arrival time ti. 

3. Justification and implementation of the wavefront model. The wave 
of advance model 2 has a deep mathematical basis. Such fronts are the salient 
feature of a wide class of phenomena that involve an exponential growth 
of some quantity (in this case, the population density) and its spread via 
random walk or related diffusion; see Fort and Pujol (2008) and references 
therein. This behavior is captured by the classical equation of Fisher (1937) 
and Kolmogorov, Petrovskii and Piskunov (1937), known as the FKPP equa- 
tion, of the form 

(2) ^ =7 A^i__j +V .(*, G VA0, 

where TV" is the population density, a function of position x and time t, 7 
is the growth rate of the population (the difference of the birth and death 
rates per unit area), K is the carrying capacity of the landscape (the maxi- 
mum sustainable population density) , and uq is the diffusivity, a measure of 
human mobility. For constant K, this equation has obvious solutions N = 
and N = K. Solutions of the initial value problem for this equation in a 
homogeneous domain, with a localized initial condition, have the form of a 
wave of advance, N(x — Ut) in one dimension, where the solution changes 
from N = ahead of the wavefront to N = K behind it. The width of the 
transition region (an internal boundary layer) is of order d = y^c/T and, 
in one dimension, the front propagates at a constant speed 

(3) U = 2^J^ 

independently of K. In two dimensions, the constant-speed propagation pro- 
vides a good approximation as soon as the radius of curvature of the front 



2 To avoid confusion, we note that the word "wave" does not refer here to any oscillatory 
behavior. 
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becomes much larger than d. For parameter values characteristic of the Ne- 
olithic expansion — for example, 7~ 0.02 year -1 and uq ~ 15 km 2 s -1 (see 
the discussions in Sections 3.2 and 3.4), providing U ~ 1 km s _1 — the front 
width is d ~ 30 km, and the model of a wavefront propagating at a con- 
stant speed is safely applicable at the continental scales of 100-1000 km. 
Even though we aim at an inference for the speed of the Neolithic expan- 
sion, we find it convenient to present our results with reference to the FKPP 
model. In particular, U is parameterized in terms of 7 and vq since these 
quantities admit a direct interpretation in terms of human behavior. Nev- 
ertheless, we emphasize that the inference that follows does not depend in 
any fundamental way on the FKPP equation, or on the assumptions of the 
demographic processes underlying it; the results we present may simply be 
considered as providing statistical constraints on the speed of a wavefront 
U, independently of any interpretation of the mechanisms behind this wave. 

3.1. The wavefront representation. A discrete version of the wavefront 
model, suitable for numerical simulations, is constructed as follows. We se- 
lect the point source of the population and then, at a small radius from 
this point, we arrange a small number of test points ("particles") on a circle 
around it, with the position of particle i denoted by Xj = (fa,Xi), where fa 
and Xi are the geographical latitude and longitude, respectively. This allows 
us to compute the local tangent and normal vectors of the front, in partic- 
ular, the local unit outward normal rij (in the direction of the advancing 
front). At each time step we move each particle along rij with the local ve- 
locity u given by (6). The wavefront thus expands outward in time; in a 
homogeneous, isotropic domain, the expansion would take the form of con- 
centric circles, but in the heterogeneous and anisotropic domain used here 
(described in detail below), the expansion is irregular and tracks the front 
obtained from the corresponding FKPP model. At the continental scale, the 
Earth's curvature is noticeable and we work on a spherical surface of the 
appropriate radius. 

The separation of the particles increases as the front expands, and an 
additional particle is introduced between any two particles if their separation 
exceeds a specified value 5. Conversely, if the front contracts, particles that 
are closer than 5/2 are removed. We take 5 to be equal to the grid spacing 
(introduced below) of 4 arc-minutes, that is, between 3 km at 60°N and 
7 km at 30°N. It is important to monitor the particle ordering along the 
front after adding or removing particles in order to maintain a continuous 
wavefront. 

When the front encounters an obstacle (e.g., a mountain range), its flanks 
move round it and then have to merge behind the obstacle. The merger is 
implemented by applying algorithms initially used to model the evolution of 
magnetic flux tubes in astrophysical simulations, and so described in more 



POPULATION DYNAMICS IN THE NEOLITHIC 



7 



detail in Baggaley et al. (2009). Briefly, we check, at each time-step, if any 
two particles (which are not neighbors) are closer than 5; if so, the particle 
ordering is switched to "short-circuit" any (almost closed) loops that have 
arisen. The front then propagates further, leaving behind closed loops which 
are then removed from the simulation since they do not affect the front 
arrival time. 

3.2. Regional variations. The front speed (3) depends on the product of 
7 and i/q , so that these parameters cannot be estimated simultaneously from 
studies of the front propagation. The growth rate is rather well constrained 
by biological and anthropological factors: Birdsell (1957) reports values in 
the range 0.029-0.035 year -1 and Steele, Adams and Sluckin (1998) suggest 
the range 0.003-0.03 year -1 . Here, for comparability with our own recent 
studies [e.g., Davison et al. (2006)], we take 7 ~ 0.02 year -1 , corresponding 
to a generation time of 30 years. This choice does not limit the applicability 
of our model. If, for example, a slightly different value of 7 is preferred, then 
the reported values of z/q need only to be scaled so that the product 71^ 
remains unchanged. 

To make the problem tractable, we adopt a specific spatial variation of the 
diffusivity and seek inference for its magnitude, just a single scalar parame- 
ter. Such assumptions unavoidably involve significant arbitrariness; however, 
a well-informed assumption is justifiable if it results in a model that fits the 
observations to a satisfactory degree. Thus, we represent the diffusivity in 
the form 



where v is the constant-dimensional magnitude and z/l is a dimensionless 
function of position whose magnitude is of order unity and whose form is 
assumed to be known. We prescribe z/l to depend solely on the topography 
(altitude above the sea level) and latitude as described below. 

Early farming does not appear to have been practical at high altitudes 
in the temperate zones considered here [e.g., sites in the Alpine foreland 
are not found on land above 1 km: Whittle (1996)]; to include this effect 
in our model, we smoothly reduce z/l to zero at around this height. Also, 
in order to allow for coastal sea travel, z/l exponentially decreases with the 
distance c?l to the nearest land. Finally, we include a linear decrease of z/q 
toward the northern latitudes. Altogether, we adopt the following form for 
the dimensionless diffusivity: 



where a is the altitude, 4> is the latitude, and we recall that a < corresponds 
to the seas. Both a and d\, are complicated functions of position reflecting 



(4) 





a > 
a < 
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the topography and coastlines of Europe. Of course, the choice of these 
particular forms are only a crude and exploratory attempt to capture the 
dependence on latitude and altitude. However, we have found our results 
not to be too sensitive to minor changes in these functions. 

The altitude data have been taken from the ETOPOl 1-minute Global 
Relief database [Geophysical Data System (2011)]. As a reasonable com- 
promise between computational efficiency and spatial accuracy, we use a 
740 x 1100 grid with a spatial resolution of 4 arc-minutes between 15° W 
and 60°E in longitude and 25°N and 75°N in latitude. The boundaries of 
the domain explored and the spatial variation of are shown in Figure 1. 

Significant regional anomalies revealed by archaeological and radiometric 
evidence are the early Neolithic Linearbandkeramik (LBK, Linear Pottery) 
Culture that propagated in 5500-4900 BC along the Danube and Rhine 
rivers at a speed of at least 4 km/year [Dolukhanov et al. (2005)] and the 
Impressed Ware culture (5000-3500 BC) that spread along the west Mediter- 
ranean coastline at an even higher speed. The latter could be as high as 10 
km/year [Zilhao (2001)]. Davison et al. (2006) included directed spread (ad- 
vection) along the major river paths and coastlines into their mathematical 
model to achieve a significant improvement in the fit to 14 C dates. 

We include advection along the Danube-Rhine river system and the coast- 
lines using world map data from Pape (2004), in the form of tangent vectors 
of the coastlines and river paths, vc and vr, respectively, given at irreg- 
ularly spaced positions. To remap the data to the regular grid introduced 
above, we use the original data weighted with exp(— dy/lh km), where dy 
is the distance between the grid point and the data point, and the chosen 
scale of 15 km contains 2-5 grid separations at latitudes 30°-60°N. 

With the unit tangent vectors for the coastlines and river paths thus 
defined, the magnitudes of the respective local front velocities, Vc and Vr, 
are subject to inference. In other words, the local front velocity is given by 

(6) u = ?7n + V, 

where V = Vc sgn(n • vc)vc + Vr sgn(n • vr)vr. Since, unlike the geograph- 
ical data, the points defining the front are not necessarily at the nodes of 
the regular grid, we used bilinear interpolation from the four closest grid 
points to calculate z^l (and thus U), vc and vr at the front positions. We 
note that the accuracy of the wavefront model was tested by comparing it 
with numerical solutions of (2) for a wide range of parameter values; the 
agreement was invariably quite satisfactory with a typical deviation of 40 
years. 

3.3. Statistical model. We model the radiocarbon dates for each object 
at each site as being a noisy version of those dates predicted by the wavefront 
model, explicitly allowing for three types of error introduced in Section 1: 
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the date tij of object j at a site i (location Xj), for j = 1, . . . , mi, i = 1, . . . , n, 
is modeled as 

tij = Ti(Q) + C(xj) + aijUij + asij, 

where Ti{6) is the deterministic wavefront arrival time, a function of 6 = 
(u, Vc,Vr) t , and cjy and are independent error terms following a stan- 
dard normal distribution. Here C( x i) allows for spatially inconsistent data, 
Oij is the standard deviation of the 14 C measurement in the laboratory and 
calibration, and a is a global error term. Between them, the £ and a terms 
allow for mismatches between our model and the data: C, allows for a certain 
level of variability between the arrival times at nearby sites, while a allows 
for an additional global variability (with no spatial correlation). As well as 
the mismatch arising when a simple large-scale model is applied to a pro- 
cess with inherent variability on smaller scales, these terms will also allow 
our model to be fairly robust to any problems in our radiocarbon data such 
as the misinterpretation of sites of secondary Neolithic occupation as being 
indicative of the wave of first arrival. We note that some authors [e.g., Buck 
et al. (1991)] choose to model the start-date of a "first arrival phase" using 
an additional parameter for each site, with the explicit expectation that all 
individual samples will fall after this date. However, our wavefront model 
attempts to capture the Neolithic transition over much larger spatial and 
temporal scales and so we choose not to incorporate such detailed temporal 
effects. 

We model the spatially smooth error process C using a zero mean Gaus- 
sian process (GP) prior [Rasmussen and Williams (2006)] with covariance 
function k^(-,-), that is, 

CO ~GP(0,fc f (-,.))■ 

We impose smoothness by taking a Gaussian kernel for the covariance func- 
tion so that the covariance between spatial errors at locations x and x' is 

k^(x,x') = a^exp{— (x — x') T (x — x')/r^}. 

The parameters of this function control the overall level of variability and 
smoothness of the process, with larger values of the length scale giving 
smoother realizations and thereby down- weighting sites with different dates 
relative to nearby sites. 

In the following analysis it will be more straightforward to work with an 
equivalent model for the summary statistics (tj,0"i) in (1), namely, 

(7) ti = Ti(0) + ((xi)+aiU}i + aEi, i = l,...,n, 

where the cjj and £i are independent error terms following a standard normal 
distribution. Thus, the inferential task is to make plausible statements about 
the wavefront model parameters 6, the global error term a and the Gaussian 
process hyperparameters and r^. 
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3.4. Bayesian inference. We adopt a Bayesian approach to inference and 
express our prior uncertainty for the unknown quantities through a density 
w(9, a, a^, r^). The speed of spread reported by Ammerman and Cavalli- 
Sforza (1971), together with the growth rate discussed in Section 3.2, sug- 
gests v ~ 12 km 2 /year. The spread reported for the LBK culture [e.g., Am- 
merman and Cavalli-Sforza (1973), Gkiasta et al. (2003)] suggests Vr, ~ 5 
km/year. The spread reported for the Impressed Ware culture [Zilhao (2001)] 
would suggest Vq ~ 10 km/year. However, this study is restricted to the 
Western Mediterranean and we believe this parameter would take a smaller 
value over the European continent such as Vc — 3 km/year. We do not ex- 
pect the wavefront model to fit the data well in all regions, so we expect a 
relatively large value of a (compared with the accuracy of individual dates) , 
a ~ 500 years. The magnitude of a^, which models the variability in dates 
between spatially nearby sites, is difficult to estimate, but the maximum ac- 
curacy of laboratory radiocarbon measurements (which will probably be the 
smallest source of variability within our model) suggests a minimum signifi- 
cant value of ~ 20 years. The mean minimum separation of radiocarbon- 
dated sites suggests ~ 1.3°. We use the values above to guide our choice of 
mode for the prior distribution. In view of the uncertainties in these values, 
we make the prior rather diffuse, with independent components. Specifically, 
we use the log-normal and inverse-gamma distributions, with 

i/~LN(3.5,l), Vb~LN(l,0.5 2 ), Vr ~ LN(2.6, 1), 

tr 2 ~IG(5,10 6 ), a c ~LN(5,1.5 2 ), r ( ~ LN(2.5, 1.5 2 ). 

Let t = (ti, . . . ,t n ) T and t{9) = (t\(9), . . . ,r n (9)) T denote the observed ar- 
rival times and those from the wavefront model (evaluated at 9). Then, using 
(7), the data model is an n-dimensional normal distribution, with 

t|0,a,a c ,r c ~N n (T(0),£), 

where £ = K^(X, X) + diag(<7 2 + <r 2 ), K^(X,X) has (i,j)th element /c^(xj,Xj), 
and X = (xi , . . . , x n ) T are the locations of the radiocarbon-dated sites. Com- 
bining both sources of information (the data/model and the prior) gives the 
joint posterior density as 

7r(0, a, a,£ , rjt) oc 7t(0, a, a^, r^)7r(t\0, a, o^, r^). 

In practice, the posterior density is analytically intractable and we there- 
fore turn to a sampling-based approach to make inferences. Markov chain 
Monte Carlo (MCMC) methods can readily be applied to this problem, but 
such schemes will typically need many thousands or millions of evaluations 
of t(9), each requiring a full simulation of the wavefront model expand- 
ing across the whole of Europe. This takes around 10 seconds on a quad- 
core CPU with a clock speed of 2.67 GHz and the code parallelized using 
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OpenMP. This computational cost prohibits the use of an MCMC scheme 
for sampling from the posterior distribution of the unknown quantities. To 
proceed, we therefore seek a faster approximation of the first arrival times 
from the wavefront model. 

4. An approximation to the wavefront model. The first arrival times 
from the wavefront model can be approximated by using both deterministic 
and stochastic methods such as cubic splines and Gaussian processes, re- 
spectively. These approximations (known as emulators) are determined by 
first setting up an experimental design consisting of a number of locations 
and parameter values, and then computing the (deterministic) first arrival 
times from the wavefront model using the locations and parameter values in 
the design. Finally, the proposed emulator is fitted to this model output. In 
line with many other authors, for example, Kennedy and O'Hagan (2001) 
and Henderson et al. (2009), we favor using Gaussian processes to emulate 
the first arrival times, as they not only fit this model output exactly and 
interpolate smoothly between design points but also quantify levels of un- 
certainty around interpolated values. Rather than build a complex emulator 
whose inputs are both site location and wavefront model parameters, we 
pragmatically build individual (parameter only) emulators for each site and 
rely on a spatial Gaussian error process to smooth the resulting spatial inho- 
mogeneities. This strategy is viable because the inference scheme described 
later in Section 5 only requires emulation at the 302 sites in the radiocar- 
bon data set. One possible drawback is that it does produce site-specific 
emulators that are not as spatially smooth as the more complex emulator. 
However, there are significant computational benefits of building separate 
emulators for each of the 302 sites, as they work in a smaller input space 
and can be computed in parallel. 

4.1. A Gaussian process emulator. The output from a single simulation 
of the wavefront model with input parameter values = (v, Vc, Vr) t gives 
the front arrival time for all sites. Consider the arrival time Tj(0) at a single 
site i. We model the emulator for the arrival time at this site using a Gaussian 
process with mean function rrii(-) and covariance function that is, 



A suitable form for the mean function can be determined by noting that the 
great-circle distance d from the source of the wavefront is approximated by 
d = ut, where u is the front speed given by (6) and r is the arrival time. 
Hence, r oc 1/u and therefore we seek 



(8) 




mi(0) = Ctifi + Oi,l—7= + a i,2T7~ + a i,37J- 




12 



A. W. BAGGALEY ET AL. 



where are coefficients to be determined, which account for the relative 
importance of diffusivity, coastal and river speeds. We determine them using 
ordinary least squares fits. Note that the mean function depends on ^Jv to 
be consistent with (3). We use a stationary Gaussian covariance function 

(0, 6') = a 2 t expj - - O'jf/rh j , 

where r^j is the length scale at site i along the fy-axis, as is widely used 
in similar applications since it leads to realizations that vary smoothly over 
the input space [Santner, Williams and Notz (2003)]. We describe how to 
make inferences about the parameters (a,, r^) of these site-specific covariance 
functions in Section A. 2. 

Suppose that p simulations of the (computationally expensive) wave- 
front model are available to us, each providing the arrival time at each 
radiocarbon-dated site. Let Tj(0) = (rj(0i), . . . , Ti(O p )) T denote the p- vector 
of arrival times resulting from the wavefront model with input values @ = 
(01, . . . , 6 P ) T , where 0j = (fj, Vc,i, VR,i) T - Prom (8) we have 

T i (e)~jv p (m i (e),iif i (e j e)) j 

where m^(O) is the mean vector with jth element mj(0,) and -fQ(0,0) is 
the variance matrix with (j,£)th element ki(Oj,0g). 

We can use the simulations to model the front arrival times corresponding 
to other values of the input parameters. Indeed, suppose we need to simulate 
front arrival times at a collection of p* new design points O* = {0\, . . . , 6**) T . 
This is straightforward as, with the local emulators, the arrival times are 
mutually independent and the distribution of the arrival time Tj(0*) at site 
i conditional on the training data Tj(O) can be determined using standard 
properties of the multivariate normal distribution as 

(9) r i (9*)|r i (e)~iv p ,( / x l (e*),y i (e*)), 

where 

(10) a*»(©*) = m 4 (e*) + Ki(e*,e)Ki(e, er^nte) - m^e)}, 

(11) Vi(e*) = Ki{@*, e*) - Ki(e*,e)Ki(e, e ; a h ^r^e, e*). 

Note that, to simplify the notation, we have dropped the dependence in these 
expressions on the hyperparameters (aj,rj), where = (^1,^2,^3)^. 

5. Inference. Before we can fit the Gaussian process emulators, we need 
to obtain the training data from the wavefront model. These data are ob- 
tained by running the emulators at a particular experimental design; see 
Section A.l for details. Given the 14 C arrival times and the training data, 
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it is possible in theory to fit the emulator and statistical model (7) jointly 
using an MCMC scheme. However, this is likely to be extremely compu- 
tationally expensive, and we therefore fit the emulator and the statistical 
model separately. This approach is advocated by Bayarri et al. (2007) and 
adopted by Henderson et al. (2009), among others. Details on how we fit 
the emulators and demonstrate their accuracy can be found in Sections A. 2 
and A. 3. 

We are now in a position to build an inference algorithm using the site- 
specific stochastic emulators t*{6) instead of the deterministic wavefront 
model Ti(9). In doing so, we fix each emulator's hyperparameters (aj,r,) 
at their posterior mean. It is possible to undertake inferences allowing for 
uncertainty in the hyperparameters. We describe methods for doing this in 
Section A. 4 and also show that allowing for such uncertainty in our analysis 
has only a very minor affect on the posterior distribution and hence on our 
inferences. 

Using the subscript e to denote the use of the emulators, the (joint) pos- 
terior density for the unknown parameters in the statistical model, replacing 
(7), is 

(12) 7r e (0,cT,a c ,r ? |t) oc 7r(0, a, o c , r^)ir e (t\0 , a, a r c ), 

where 

7r e (t\e,a,a r c ) cc |X(0)p 1/2 exp{-±(t - n(e)) T E(G)- 1 (t - p(0))}, 

/x(0) has ith element /J«(0) and, to allow for emulator uncertainty, we re- 
define E(0) = K C (X,X) + diag{Vi(0) + of + cr 2 }, where X = (xi, . . . ,x n ) T . 
Sampling from the posterior distribution (12) can be achieved by using a 
Metropolis-Hastings scheme similar to that described in Section A. 2. Specif- 
ically, we use correlated normal random walks to propose updates on a log- 
scale in separate blocks for 6, a and the parameters (a^,r^) that govern 
spatial smoothness. The covariance structure of the innovations was chosen 
using a series of short trial runs. 

6. Results. We found that running the MCMC scheme for 1M iterations 
after a burn-in of 100K iterations and thinning by taking every 100th iterate 
gave a sample of J = 10K (effectively uncorrelated) values from the posterior 
distribution. Figure 2 shows the marginal and joint posterior densities of the 
key parameters of the wavefront model (the diffusivity v, the coastal velocity 
Vc and the river-path velocity Vr) and also the parameters of the statistical 
model (the global error a, and the Gaussian Process hyperparameters 
and 77;). The reduction in uncertainty from the prior to the posterior, which 
is considerable for most parameters, shows that the data have clearly been 
informative. 
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Fig. 2. Marginal (upper row, from left to right) and joint (middle row) posterior densi- 
ties for the wavefront model parameters v , Vc and Vr, and marginal posterior densities 
(bottom row, left to right) for the residual uncertainty parameter o and the spatial GP 
hyperparameters and rj. The prior densities are shown in dashed red. 



The posterior distributions of the wavefront model are in reasonable agree- 
ment with the range of values suggested by previous studies. For exam- 
ple, the mode of the marginal posterior for the diffusivity u is around 20 
km 2 /year, which corresponds to a global rate of spread of U = 2^/tz7~ 1.25 
km/year. The modal values for the amplitudes of the advective velocities, 
Vc = 0.3 km/year and Vr = 0.2 km/year, are rather lower than the val- 
ues suggested in the archaeological literature which motivated their study 
[Zilhao (2001), Dolukhanov et al. (2005)]. Regarding the coastal speed Vc, 
this could be explained by an accelerated spread only being required in the 
west Mediterranean, whereas our model applies it to the whole European 
coastline; an alternative model, allowing regional variations in coastal ef- 
fects, may produce different results. The relatively small value for Vr is 
perhaps more surprising, given the prevailing archaeological picture of rel- 
atively rapid spread of the LBK culture in the Danube and Rhine basins. 
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It may be that the earliest dates associated with the LBK culture do not 
correspond particularly well to spread as a continuous wave and that an 
alternative model might also be preferable here. After taking account of the 
overall level of fit of our current model, the dates within this region are 
satisfactorily explained with a relatively weak advective enhancement of the 
background diffusive spread. It is also plausible that the slightly enhanced 
value of v (compared with many of the studies cited in Section 1) results in 
the relatively low values for Vc and Vr. 

The marginal posterior distribution for the global error parameter a is 
centred on a value of around 575 years (see Figure 2). This value quantifies 
the mismatch between our mathematical model of the spread and the local 
variations present in the actual spread. Thus, our inference suggests that a 
simple large-scale wave of advance, while remaining a good model on the 
continental scale, should only be considered a good model on time scales 
of order 600 years (and thus length scales of order 600 km) or greater; on 
shorter time scales, significant local variations should be expected. Mismatch 
between model and data is also allowed in our spatial term, £; the marginal 
posterior distributions for its amplitude (a^) and length scale (Vf) are rela- 
tively tightly peaked around values of order 30 years and 10°, respectively. 
The latter value suggests that correlations between dates at nearby sites 
are effectively significant within distances of order 750 km. (Note that this 
is approximately the length scale on which a suggests that the variations 
should be expected within the large scale model, so that the two results are 
consistent.) However, the amplitude is relatively small compared with a, 
suggesting that "within region" variation is not particularly significant, in 
comparison with larger scale deviations from the model. 

We can examine the operation of the spatial Gaussian process in more 
detail by looking at its distribution at various sites. Now 

(13) C|t, 0, a, a r c ~ N n (pjX), V*(X)), 

where 

H.(X) = K C (X, X) T {K C (X, X) + E*} -1 {t - »(0)}, 

V*(X)=K ( (X,X) - K ( (X,X) T {K C (X,X) + E,}~ 1 K C (X,X), 

and so posterior realizations from the spatial process can be obtained us- 
ing the MCMC output (0^\a^\a^\r^), j = 1, . . . , J, by simulating from 

^\t,0^\a^\a^\r^ , j = 1, . . . , J. For example, the posterior spatial pro- 
cess at the UK site Monamore (5.1°W, 55.5°N) has mean 21.7 years and 
standard deviation 33.0 years, and this small mean reflects both the good fit 
of the wavefront model at this site and at its neighboring sites; see Figure 1. 
In contrast, the French site Greifensee (8.7°W, 47.7°N) has mean -726 years 
and standard deviation 453 years, and this large mean highlights that the 
data from this site is inconsistent with its neighbors. 



16 



A. W. BAGGALEY ET AL. 




Fig. 3. 14 C dates (red squares, displaced vertically for presentation purposes) and 
their predictive densities (solid) for four randomly chosen sites: (a) Monamore 
(mi = l,ai = 120), (b) Vrsnik-Tarinci (mi — 22, cr, = 13), (c) Galabnik (nii = 10, en = 16) 
and (d) Greifensee (mi = 1,0^ = 130/ 



6.1. Model fit. We use predictive simulations to assess the validity of 
the statistical model and thereby the underlying model of the propagating 
wavefront. The (joint) posterior predictive distribution of the arrival times 
tprcd at the n = 302 radiocarbon sites can be determined in a similar way to 
that for the (posterior) spatial process. First we have 

T'pred 

where £*(0) = diag{T^(0) + of + a 2 } and £ = (C(xi), • • • , C(*n)) T is the spa- 
tial Gaussian process evaluated at the sites. We can integrate over £ in closed 
form using (13) to obtain 

t prcd 1 1 , 6 , a, a c , r c ~ N n (n(G) + ^ (X) , £* (0) + K (X) ) 

and thereby determine the predictive arrival time distribution by simulat- 
ing realizations from t^ d |t, a^, r£ , j = 1, . . . , J. Figure 3 shows 
the radiocarbon dates for four randomly chosen sites together with their 
posterior predictive distributions determined using this simulation tech- 
nique. These sites are at Monamore (5.1°W, 55.5°N), Vrsnik-Tarinci (22.0°E, 
41.7°N), Galabnik (23.1°E, 42.4°N) and Greifensee (8.7°W, 47.7°N). The 
14 C dates cluster very close to the mode of the posterior densities in three 
cases and the discrepant site (d) is the one discussed in the previous section 
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Fig. 4. TTie mean (left) and standard deviation (right) of the posterior predictive distri- 
bution of the arrival time. 

whose data is not consistent with its neighboring sites. Very encouragingly, 
we found similar clustering of site dates around their predictive mode for 
around 90% of the 302 sites. We conclude that the wavefront model (or 
more properly, the emulators) provides an adequate description of the data. 
Figure 4 illustrates how the wavefront moves across Europe via a map of the 
means of the posterior predictive distributions (Figure 3) generated by inter- 
polating the mean arrival time at each site onto a regular mesh, using cubic 
interpolation with the Matlab routine griddata. In the same figure we also 
display the standard deviation of the predictive distributions interpolated 
in the same manner, highlighting areas where there is large uncertainty in 
the timing of the wavefront arrival from our model. 

7. Discussion. Our main goal in this paper has been to develop generic 
statistical tools rather than to propose a more complex model for the spread 
of the Neolithic than, for example, Ackland et al. (2007), Davison et al. 
(2009) and Isern and Fort (2010). Our approach to inferring parameters of 
the Neolithic expansion is based on a direct application of the wavefront 
model. To alleviate the computational demands of the MCMC scheme, we 
constructed Gaussian process emulators for the arrival time of the wavefront 
at each of the 302 radiocarbon-dated sites in our data set, with the input 
parameters obtained from the wavefront model. This approach has been 
shown to perform well and is likely to do so for other spatio-temporal models 
that have complex space-time dependencies and are slow to simulate. 

This paper represents the first attempt of statistical inference for quantita- 
tive models of the Neolithic dispersal based on modern statistical methods, 
providing an opportunity for a more detailed and reliable analysis of the 
data than ever before. We have improved upon earlier model fitting for the 
Neolithic expansion by avoiding the introduction of a maximum precision 
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of the Early Neolithic C age determinations at 100-200 years, obtained 
as the minimum empirical standard deviation of 14 C dates in large homo- 
geneous data sets for key Neolithic sites [Dolukhanov and Shukurov (2004), 
Dolukhanov et al. (2005), Davison et al. (2009)]. In fact, results presented 
here remain unaffected when such a maximum precision is introduced. This 
was not the case with the earlier analyses, and the robustness of the inference 
procedure suggested here is encouraging. 

We believe that we have demonstrated the usefulness of our method as 
a generic tool and that adopting a statistically rigorous approach consti- 
tutes a significant advance on previous "wave of advance" models of the 
Neolithic. The posterior distribution of the parameters of our wavefront 
model (^, Vc,Va) provides clearly bounded constraints on the human pro- 
cesses underlying the spread of farming. Also, the posterior distribution of 
the parameters of our statistical model (a,a^,r^) allows a quantitative as- 
sessment of the utility of our wavefront model (including the limits of appli- 
cability of such a large-scale model to smaller regions) . The latter point has 
important implications for the archaeological interpretation of the Neolithic 
transition in Europe and should contribute to the ongoing debate about the 
characterization of the spread as a continuous phenomenon. For example, 
Rowley-Conwy (2011) argues that the expansion of farming was character- 
ized by sporadic events, punctuated by long periods without outward expan- 
sion; Rowley-Conwy includes the spread of the LBK and Impressed Ware 
cultures as rapid events of this type. In light of the difficulties encountered 
here in modeling these spreads, noted in the previous section, it would be 
of interest to consider alternative models, which are perhaps better able to 
model the expansion within the relevant regions. 

In terms of future work, the wavefront model can be improved in several 
respects. Most importantly, a more detailed description of regional variations 
appears to be necessary. In particular, 14 C data suggest that the accelerated 
spread along the coastlines is most pronounced in the western Mediterranean 
and around the Iberian Peninsula, but not on the North Sea and Black Sea 
coasts. The approach presented here can be readily generalized to include 
such a variation and to provide a useful estimate or an upper limit on the 
coastal propagation speed in Northern Europe. Also, our understanding of 
the Neolithic spread can be enriched by adding other parameters to our 
model, such as the starting time and position of the source of the Neolithic 
expansion, and even additional sources as suggested by Davison et al. (2009), 
based on the 14 C evidence and corroborated by Motuzaite-Matuzeviciute, 
Hunt and Jones (2009) using archaeobotanical evidence. Another major im- 
provement would be to use the published 14 C dates directly for inference, 
rather than combining them at each site as in (1). Model fits to the original 
radiocarbon data without preliminary selection, unavoidably involving sub- 
jective and possibly poorly justified criteria, have never been attempted in 
the studies of the Neolithic [Hazelwood and Steele (2004), Steele (2010)]. 
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A perennial question in the ongoing debate about the nature of the Ne- 
olithic is the mechanism of the spread of the agropastoral lifestyle. The 
wavefront model used here applies equally well to demic spread, involv- 
ing migration of people, and to cultural diffusion, that is, the transition 
to farming due to indigenous adaptation and knowledge transfer [Zvclcbil, 
Domam-ska and Dennell (1998)]. However, cultural diffusion may be mod- 
eled more explicitly, in terms of multiple, mutually-interacting populations; 
see, for example, Aoki, Shida and Shigesada (1996), Ackland et al. (2007) 
and Patterson et al. (2010). We intend to explore the fit of some of these 
other models and to determine how plausible these different models are as 
explanations of the radiocarbon data. 

APPENDIX 

A.l. Experimental design. In order to fit the emulators, we first need 
to obtain the training data by running the wavefront model at a suitable 
choice of 0-values. Computing resources available to us allowed a p = 200- 
point design. Although using a regular lattice design is appealingly simple, it 
is not particularly efficient. Instead, we adopt a more commonly used design 
for fitting Gaussian processes, namely, the Latin Hypercube Design (LHD) 
[McKay, Beckman and Conover (1979)]. Designs of this class are space-filling 
in that they distribute points within a hypercube in parameter space more 
efficiently than a lattice design. Our 200-point LHD was constructed using 
the maximin option in the MATLAB routine lhsdesign. We set the lower 
bounds of the hypercube to be the origin and used the upper one percentiles 
of the prior distribution as its upper bounds. The routine then simulated 
5000 LHDs and took the design which maximized the minimum distance 
between design points. We then used a sequential strategy in which we re- 
peatedly ran the inference algorithm (described in the following section), 
used the results to determine a conservative estimate of a hypercube con- 
taining all points in the MCMC output (and therefore plausibly containing 
all of the posterior density), and then generated another LHD. The final 
LHD used for inferences on (v, Vc,Vr) t in this paper is contained within 
the hypercube (0,120) x (0,3) x (0,2). 

A. 2. Fitting the emulator. At site i, this essentially requires the determi- 
nation of the posterior distribution of the emulator's hyperparameters using 
the training data Tj(0) from the wavefront model evaluated at the emula- 
tor's design points. We adopt fairly diffuse independent log-normal priors for 
the GP hyperparameters, with en ~ LN(6.3, 0.5) and r^j ~LN(9, 3). Then 
their (joint) posterior density is 

ir(ai,ri\Ti) oc ^(a^vr^i)^^^)^^^)^^!^, r^), 
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Fig. 5. Solid (black) curves: the marginal posterior densities of the emulator's hyperpa- 
rameters ai and ri (from left to right) for the Achilleion site. Dashed (red) curves: the 
corresponding prior probability densities. 

where 

7r(Ti|ai,rt) oc |Kj(ai,rj)|~ 1/2 exp{-i(ri - vcn) T , ,K i (a i ,r i )~ 1 (Ti - mi)}. 

For notational simplicity, we have dropped the explicit dependence of the 
training data Tj , the mean m, and the covariance matrix Ki on the training 
points 0. It is fairly straightforward to implement an MCMC scheme to 
simulate realizations from this posterior distribution. We used a Metropolis- 
Hastings scheme in which proposals for each hyperparameter were made on 
a log-scale using (univariate) normal random walks. A suitable magnitude 
of the innovations (their variance) was adopted after a series of short trial 
runs. For each site-specific emulator, the MCMC scheme was run for 500K 
iterations after a burn-in of 50K iterations and the output thinned by taking 
every 50th iterate, giving a total of 10K (almost uncorrelated) iterates to 
be used for inference. Figure 5 shows the marginal posterior densities of the 
emulator hyperparameters for site 1, the Achilleion site, located at latitude 
39.2°N, longitude 22.38°E. Note that these posterior densities have relatively 
small variability and are robust to the parameter choice for the prior of the 
hyperparameters, and both these features are typical of all the sites in our 
data set. 

A. 3. Emulator validation. The accuracy of our fitted emulators as an 
approximation to the wavefront model can be assessed by using a variety 
of graphical and numerical methods. Our assessment was based on the fit 
at a set of parameter values that were not used to fit the emulator: we 
used a p* = 100-point Latin hypercube design 0* = (0*, . . . , 0*») T . First, for 
each site i, we obtained the arrival times r« = (xj(0*), . . . , Tj(0*»)) T from 
the propagating front model. We then obtained realizations of the emulator 
at these design points which properly accounted for the uncertainty in the 
hyperparameters by using (9) and the realizations from the MCMC fit of the 
emulator. Plots of the emulator's mean and 95% credible region showed that 
the emulators provided a reasonably accurate fit to the wavefront output. 
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Fig. 6. Diagnostic plots assessing the fit of the emulator to the wavefront model. Left 
panel: the Mahalanobis distance MDi for each site, together with its distribution's upper 
95% point (dashed red line at MDi = 11.44J. Right panel: the spatial distribution of sites 
colored according to the Mahalanobis distance. 

We note that these plots were very similar to those determined by fixing the 
hyperparameters at their posterior means. 

We also assessed the fit of our emulators by using a numerical statistic 
taken from Bastos and O'Hagan (2009) which compares the emulator output 
with that from the wavefront model allowing for the site-specific accuracy 
and correlation between design points. This statistic is the site-specific Ma- 
halanobis distance MDj, i = 1, . . . , n, where 

MDf = ( Ti - ^(e*)) T v-(e*)- 1 (r i - Mi (e*)), 

which has a scaled .F-distribution in the case of the GP emulator with fixed 
hyperparameters, with MD- ~ p*(p — 5) F p * >p _3 / (j> — 3). In our calculations 
we have fixed the hyperparameters at their posterior means. Figure 6 shows 
the Mahalanobis distance at each site, together with the upper 95% point 
of its distribution, and the spatial distribution of the Mahalanobis distance. 
The figure confirms that the emulators provide a reasonable fit throughout 
the design space and that there is no obvious spatial dependence in the fit 
of the emulators. We note that these diagnostics gave similar conclusions for 
all values in the posterior sample of (aj,r,) and for different LHDs without 
any systematic site-specific biases. 

Finally, we have assessed the distributional fit of the emulators by cal- 
culating the probability integral transform (PIT) values [Gneiting, Balab- 
daoui and Raftery (2007)] for each site at each of the 100 points in the 
design. Again, we fix the hyperparameters at their posterior means. A well 
fitting emulator should produce values from the standard uniform distri- 
bution. A histogram of the PIT values over all sites is shown in Figure 7. 
The figure also shows a summary of how closely the PIT values from each 
site resemble the standard uniform distribution. Here we use a Pearson x 2 
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Fig. 7. Analysis of probability integral transform (PIT) values, calculated for each site 
from 100 emulated arrival times t* . Left panel: a histogram of all PIT values from each of 
the 302 sites. Right panel: the \ 2 statistic, computed at each site, assessing the goodness 
of fit of the PIT values to a uniform distribution. The 95% point of the \g distribution is 
plotted as a dashed line. 



statistic after binning the data into ten equal sized bins and so also show 
the 95% point of the Xg distribution. This figure confirms that the emulator 
fits are satisfactory. 

A. 4. Allowing for hyperparameter uncertainty. In Section 5 we describe 
how to obtain realizations from the posterior distribution (12) when fixing 
the emulator hyperparameters at their posterior mean. Here we extend the 
method to take account of hyperparameter uncertainty by using their re- 
alizations {(dp , r[ = 1, . . . ,n e ,i = 1, . . . , n} obtained when fitting the 
emulators. 

We first describe an MCMC scheme in which the emulators themselves 
are marginalized over hyperparameter uncertainty. In this case, the marginal 
distribution of t*(0) can be well approximated by an equally weighted nor- 
mal mixture, namely, 

■ , Tie 

where /Xj(-;a, r) and Vf(-;o,r) are given by (10) and (11). Therefore, the 
(emulator) likelihood in (12) can be written as 

He -. 

7r e (t|0,cT,a c ,r ? ) = V-</> n (t| / x(0;aW,r^),S(0;aW,rW)), 
j'=i 

where </> n (-|/z, S) is the density of an n-dimensional normal distribution with 
mean fi and covariance matrix S, /z(0;a, r) has ith element fii(0;ai,ri) and, 
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Fig. 8. Marginal posterior densities for the wavefront model parameters (y, Vc,Vr) and 
for the statistical model parameters (u, a^,r^) allowing for uncertainty in the hyperparam- 
eters (black curves) and for fixed hyperparameters (dashed red curves). 



to allow for the emulator uncertainty, we redefine E(0;a, r) = K^(X,X) + 
di&g{Vi(0;ai,r,i) + af + a 2 }, where X = (xi, . . . , x n ) T . Sampling from the 
posterior distribution when using this (emulator) likelihood can be achieved 
by using the same methods as described in Section 5. However, a major 
drawback of using this scheme is the time taken to calculate the (emu- 
lator) likelihood at each iteration. Therefore, we prefer to use an MCMC 
scheme which uses a joint update for the model parameters and the hy- 
perparameters and has (12) as a marginal distribution. Specifically, in the 
joint update, the value proposed for the hyperparameters is sampled from 
the posterior realizations described above. Although this scheme has a much 
larger state space, each MCMC iteration is considerably faster than in the 
marginal scheme and we have found it to have good convergence proper- 
ties. 

We ran the MCMC scheme for 1M iterations (after a burn-in of 100K 
iterations), and then thinned the output taking every 100th iterate to leave 
an effectively uncorrelated sample of 10K draws from the posterior densities. 
Figure 8 shows the marginal posterior distributions allowing for emulator 
hyperparameter uncertainty and for fixed hyperparameters. We note that 
it is difficult to distinguish between these cases. Similar conclusions were 
found (when comparing fixed with uncertain hyperparameters) for the site- 
specific posterior distribution of the spatial process £|t and for the predictive 
distributions i pre d|t- 
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